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Abstract 

Lattice parameter, electronic structure, mechanical and thermodynamic properties of ThN are 
systematically studied using the projector-augmented-wave method and the generalized gradient 
approximation based on the density functional theory. The calculated electronic structure indicates 
the important contributions of Th 6d and 5/ states to the Fermi-level electron occupation. Through 
Bader analysis it is found that the effective valencies in ThN can be represented as Th +1 ' 82 N -1 ' 82 . 
Elastic constant calculations shows that ThN is mechanically stable and elastically anisotropic. 
Furthermore, the melting curve of ThN is presented up to 120 GPa. Based on the phonon dis- 
persion data, our calculated specific heat capacities including both lattice and conduction-electron 
contributions agree well with experimental results in a wide range of temperature. 
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I. INTRODUCTION 



Actinide nitrides have been extensively studied in experiments in connection with their 

fl 

potential applications in the Generation-IV reactors These reactors raise a number of 
concerns surrounding the issue of nuclear power. The effective utilization of nuclear power 
will require continued improvements in nuclear technology, particularly related to safety 
and efficiency Nowadays, except for the oxide based fuels, the nitride fuels also participate 
in the competition to become the alternative materials for their superior thermophysical 
properties, such as high melting point, high thermal conductivity, and high density, as well 
as the good compatibility with the coolant (liquid sodium) [2j-|5j. Since the high density of 
nitride fuels can bring out more excess neutrons, therefore, they have a higher potential to 
transmute the long lived fission products. As for reprocessing feasibility, actinide nitrides 
also appear to be a compromise between oxide and metal fuels. For the sake of better 
understanding of the behavior of these materials under irradiation, their accurate electronic 
structure description by first-principles methods is necessary. Actinides form an isostructural 
series of mononitrides with a simple rock-salt type structure and a complete solid solubility 
in the whole composition range. In virtue of the prospective use of actinide mononitrides 
as advanced fuel materials, it is of crucial importance to know their thermal properties for 
modeling the fuel behavior at elevated temperatures. The thermodynamic properties such 
as standard enthalpies of formation and heat capacities are essential to predict the phase 
stability including the melting points. 

Despite the abundant research on actinide mononitride, however, for ThN compound, 
relatively little is known regarding its chemical bonding, mechanical properties, high pres- 
sure melting points, and phonon dispersion. In the early theoretical studies, only the lattice 
parameters and bulk modulus of the actinide nitride series have been calculated by utiliz- 
ing .inear muffin-tin orbita! (LMTO) method g Although the electronic properties and 
chemical bonding in ThN have been recently calculated by Shein et al. [7| through the full- 
potential linear-augmented-plane-wave (FLAPW) method, the study of the bonding nature 
of Th-N ionic/covalent character is still lacking. These facts, as a consequence, inhibit the 
fundamental understanding of thorium momonitride from basic point of view. Motivated 
by these observations, in this paper, we present a first-principles density functional the- 
ory (DFT) study by calculating the structural, electronic, mechanical, and thermodynamic 
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properties of ThN. As will be shown, a preliminary study of some bulk properties of ThN 
indicates that the conventional generalized gradient approximation (GGA) for the exchange- 
correlation potential in DFT can give satisfactory results when compared to experimental 
data and thus does not require other treatments beyond LDA/GGA such as LDA/GGA+ U, 
which are indispensable in the uranium and plutonium compounds for stronger localization 
and correlation of their 5/ electrons. 

After optimizing the ground state structure of ThN, we performed the calculation of 
density of states (DOS) and Bader analysis 8] for ThN. The results indicate that the Th 5/ 
states are involved in the formation of Th-N and Th-Th interatomic bonds and about 1.82 
electrons transfer from each Th atom to its surrounding N atoms. Through the mechanical 
analysis, we find that NaCl-style ThN is mechanically stable and elastically anisotropic. 
In order to predict the melting curve, we calculated the elastic constants, bulk modulus 
and shear modulus from ambient pressure to 117 GPa. Since we have known that ThN 
melts congruently at 2790±30° under a nitrogen pressure somewhat less than 1 atm 
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thus the ThN can be ranked as refractory materials. Utilizing the Lindemann melting 



criterion [lOj], we got the melting curve versus pressures. The melting points were enhanced 
by about 2100 K from ambient pressure to 117 GPa. Our predicted results indicate that 
ThN is able to withstand temperatures above 5100 K without chemical change and physical 
destruction under high pressures. The calculated phonon dispersion confirms the dynamic 
stability for ThN. Based on our phonon dispersion data, the lattice vibration energy, thermal 
expansion, and specific heat are obtained by using the quasiharmonic approximation (QHA). 
Our calculated special heat, including both lattice and conduction electron contributions, 
agrees well with experimental results at T<1500 K domain. 

The rest of this paper is organized as follows. The first-principles computational details 
are briefly introduced in Sec. II. The calculated results are presented and discussed in Sec. 
III. Finally, a summary of this work is given in Sec. IV. 

II. COMPUTATIONAL METHODS 

The first-principles total energy calculations were carried out using the Vienna ab initio 



simulations package (VASP) [U[ with the projected-augmented-wave (PAW) pseudopoten- 



tials [12] and plane waves. The exchange and correlation effects were described within GGA 
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TABLE I: Strains used to calculate the elastic constants of NaCl-type ThN. 



Strains Parameters (unlisted: £ij=0) y 9 affi*^ \s=o 

1 e xl =8 Cu - P 

2 8^=822=5 2(C n +C 12 -P) 

3 £13=^31=^ 4C 44 - 2P 



13l | . The thorium 6s 2 7s 2 6p 6 6d l 5f 1 and nitrogen 2s 2 2p 3 electrons were treated as valence 



electrons. The electron wave function was expanded in plane waves up to a cutoff energy 
of 500 eV. We have performed numerous convergence tests on determining the influence of 
the fc-point mesh on the total energy. The Monkhorst-Pack [14J llxllxll mesh was used 
in Brillouin zone (BZ) integration, which turns to be sufficient to get results converged to 
less than l.OxlO -4 eV per atom. The corresponding electronic density of states (DOS) was 
obtained with 19x19x19 /c-point mesh. 

It is known that the elastic constants are defined by means of a Taylor expansion of the 
total energy, E( V,S), for the system with respect to a small strain 5 on the equilibrium cell 
according to the following law [15 ]: 



E(V,6) = E{V ,0) + V 
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(1) 



where E(V , 0) and V are the total energy and volume of the equilibrium cell without 
strains, respectively, Tj is an element in the stress tensor, and is a factor presenting to 
take care of the Voigt index [161 ]. 

To calculate the elastic constants of rock-salt structure ThN, we applied three independent 
strains. The parametrizations that we used for these strains are given in Table [B We 
calculated the total energy of each strain for a number of small values of S. These energies 
were then fitted to a polynomial in 5 and the curvatures of the energy versus 5 curve were 
obtained for using in Eq. (1). The elastic constants are closely related to many physical 
properties of solids, such as the Debye temperature, specific heat, melting temperature, and 
Griineisen parameter, etc. At low temperatures, the vibrational excitations arise solely from 
acoustic vibrations. Hence, the Debye temperature calculated from elastic constants is the 
same as that determined from specific heat measurements at low temperature. The relation 
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between the Debye temperature (9n) arid the average elastic wave velocity (v m ) is 

where h and ks are Planck and Boltzmann constants, respectively, n is the number of 
atoms in the molecule and Vt is molecular volume. The average elastic wave velocity in the 
polycrystalline materials is approximately given by 

-1/3 



(3) 



-(— + — 

3 vf vf j 

where Vt—yGjp (p is the density) and vi=\J (3B + 4G)/3p are the transverse and longitu- 
dinal elastic wave velocities, respectively. 

After obtaining the Debye temperature at various pressures, furthermore, we performed 



the melting curve calculation 
Lindemann melting criterion 



rom Debye temperatures. The calculation was based on the 



lo| . This model is based on the harmonic approximation, 
predicting that melting will occur when the ratio of the root mean square (rms) atomic 
displacement to the mean interatomic distance reaches a certain value (generally about 
1/8). It can be expressed as follows 

T m = CV 2 / 3 9 D , (4) 
where T m is the melting point, C is a constant, and V is atomic volume. 

III. RESULTS AND DISCUSSIONS 
A. Atomic and electronic structures 

In order to examine the possibility of magnetism in ThN, the spin-polarized calculations 
by assuming initial ferromagnetic state were carried out. We found that the ground state for 
ThN is non-magnetic, without any localized atomic magnetic moments, the result of which 
is in agreement with observed [3j paramagnetism in ThN. In this paper, the theoretical 
equilibrium volume V , bulk modulus B and the pressure derivative B' are obtained by 



fitting the third-order Brich-Murnaghan equation of state (BMEOS) 17] . Our calculated 



lattice parameter for the cubic unit cell of ThN is ao=5.179 A, which is in good accordance 



with the experimental data of 5.16 A 18]. The bulk modulus B and its pressure derivative B' 




FIG. 1: (Color online) Total and site-projected orbital-resolved electronic densities of states for 
ThN at equilibrium. The Fermi energy is set at zero 

obtained by fitting the BMEOS are 176.1 GPa and 3.9, respectively, which are also consistent 
with the corresponding experimental values of 175 GPa and 4.0 [19 1. The calculated B for 
ThN is much larger than that of metallic a-Th (about 60-62 GPa 201]). i.e., a pronounced 
increase of structural rigidity from metal to nitride due to the direct Th-N bonding formation, 
which occurs as well as for other NaCl-type metal mononitrides. 

The total DOS and site-projected orbit-resolved DOS (PDOS) of ThN are displayed in 
Fig. 1. Evidently, the conduction band is strongly marked by Th 5/ states, with a little 
bit degrees of Th 6d and N 2p states. The valence band ranging from —5.5 eV to —1 eV 
is of mixed Th d/f and N p character. Near the Fermi level are composed by comparable 
contributions of Th 6d and 5/ states, whereas the contributions of N 2p states become much 
smaller. Due to the observable d and / states near the Fermi energy, ThN exhibits a clear 
metallic behavior. 

In order to further analyze the chemical bonding nature, we display in Fig. [2] the charge 
distribution in ThN (100) plane and list in Table [Til the Bader effective charges. We find 
the near-spherical distribution of electron density around thorium and nitrogen atoms, with 
rather small density value between them. This is typical for crystals with the rock-salt struc- 
ture having the ionic bonding due to the charge transfer from metal to non-metal atoms. 
Since Bader analysis is an effective tool for studying the topology of the electron density 
and suitable for discussing the ionic/covalent character of a compound, thus we also per- 
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FIG. 2: (Color online) Valence charge density of ThN in (100) plane. 
TABLE II: Calculated effective atomic charge and volumes according to Bader partitioning of AnN 



(An=Th, U, Pu). 



Compound 


Qb(A) 


Qb(N) 


V B (A) 


V B (N) 




(e) 


(e) 


(A3) 


(A3) 


ThN 


+1.82 


-1.82 


20.924 


13.848 


UN 


+1.71 


-1.71 


18.301 


11.581 


PuN 


+1.59 


-1.59 


17.713 


11.757 



formed the Bader effective charge calculation for ThN. For this we adopted 300x300x300 
charge density grids, then the spacing between adjacent grid points is 0.0172 A. The calcu- 
lated effective atomic valance charges and volumes are listed in Table [Til together with the 
UN and PuN results for comparison. Our results show that about 1.82 electrons transfer 
from each Th atom to its neighboring N atoms. The effective valency in ThN then can be 
represented as Th +L82 N" L82 , while UN is U +Ln N" Ln 21] and PuN is Pu +L59 N~ L59 , in- 
dicating that the ionicity in ThN is somewhat stronger than that in UN and PuN, as shown 
in Table II. 



B. Mechanical properties 



For the cubic structure, there are three independent elastic constants Cn, C12, and C44. 
Under pressure P, the relation between the elastic constants Cy and the bulk modulus B 




FIG. 3: (Color online) (a) Evolution with pressure of the ThN elastic constants, bulk modulus and 
shear modulus; (b) Evolution with pressure of the transverse elastic wave velocity vt, longitudinal 
wave velocity v\, and the average wave velocity v m . 



can be expressed as follows |22| : 



C u + 2C 12 = 3B — P. 



(5) 



At zero pressure, the well-known expressions of Voigt and Reuss 16|, |23|, |24| bulk moduli 
can be re-derived. Under the Voigt approximation, the effective shear modulus Gy for 
cubic phase can be expressed as Gy = Cll ~ Cl 5 2+3C ' 44 , while with Reuss approximation the 



shear modulus can be expressed as, Gr = 4c^+3(clf-^Ci2) ' —I P rove d that the Voigt 

and Reuss equations represent upper and lower limits of the true poly crystalline constants, 
and recommended that the shear modulus G is an arithmetic average of Voigt and Reuss 
approximations, i.e., G = |(Gr + Gy). From that, the Young's modulus E and Poisson's 
ratio v can be given by E= £b+g anc ^ v= 2(3B+g) • Using the above functions, the calculated 
bulk modulus B, the pressure derivative of bulk modulus B\ shear modulus G, Young's 
modulus E, and Poisson's ratio v of ThN are given in Table IIII1 For comparison, the 
experimental data from Ref. 19) and the theoretical FLAPW-GGA results in Ref. 



arc 



also presented. As can be seen from Table IHIl our calculated C\i and C44 are in agreement 
with the FLAPW-GGA results, while Cn=334.8 GPa is smaller than the corresponding data 
396.6 GPa. At ground state, the bulk modulus B derived from elastic constants is 179.4 GPa, 
which is well consistent with that obtained by BMEOS fitting. This value is very close to the 
experimental data of 175 GPa [19] but smaller than the FLAPW-GGA result of 199.9 GPa. 
The pressure derivative of the bulk modulus B' is 3.9, which is also in good accordance 



S 



TABLE III: Calculated elastic constants, elastic moduli, pressure derivative of the bulk modulus 
B' , Poisson's ratio v, and anisotropic factor A for ThN at GPa. Except B , v and A, all other 
values are in units of GPa. 



Method 


C u 


Cl2 


C44 


B 


B' 


Gy 


E 


V 


A 


GGA 


334.8 


101.7 


75.0 


179.4 


3.9 


89.5 


230.3 


0.286 


0.64 


FLAPW-GGA a 


396.6 


101.5 


79.9 


199.9 




102.4 


262.5 


0.281 


0.54 



Expt. 6 175 4.0 



26] 



19] 



with the experimental data of 4.0. The calculated Poisson's ratio i/=0.286 is consistent 
with the FLAPW-GGA result. With a cubic system, the elastic anisotropy factor is given 
by A— g ^l 4 g ■ Since that micro-cracks in materials can be easily induced by significant 
elastic anisotropy, so it is important to evaluate anisotropic factors in understanding their 
mechanical durability. Our present value of A is equal to 0.64, indicating that ThN is 
elastically anisotropic. As for isotropic crystals the factor A is equal to 1.0, while with any 
value larger or smaller than 1.10 indicating elastic anisotropy. 

Besides, Figure [3]^a) shows the pressure dependence of the elastic constants, bulk modulus 
B, and shear modulus G. One can clearly see that the elastic constants all linearly increase 
with pressure. The value of Cu, Ci 2 and C 44 are enhanced by 705 GPa, 27 GPa and 51 
GPa, respectively, from GPa to 117 GPa. These elastic constants satisfy the generalized 
elastic stability criteria for cubic crystals under pressure, i.e., C\\ > |C*i2|, C11 + 2 1 12 1 
> 0, C 44 > 0, where C 'a=Cu—P (i=l, 4) and C\i = C12+P. The dependence of bulk 
modulus B and shear modulus G on pressure also linearly increase with pressure. These 
pressure dependence of elastic constants will be used in calculating the Debye temperature, 
and further predicting the melting temperatures at elevated pressures. 



C. Melting curve of ThN 

The melting curves of materials have great scientific and technological interest. To un- 
derstand the transition between solid and liquid phases of ThN at high pressures, we will use 
the Lindemann criterion [Eq. (|4])] to calculate the melting curve. The Lindemann criterion 
will be used here as a single-parameter model and the free constant C can be calculated 
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FIG. 4: (Color online) The melting temperature of ThN as a function of pressure. 

from a single point. In order to perform the evolution of T m , we should know the melting 
temperature T m at P=0 GPa and the evolution of Op. On the one hand, we have known 
that ThN melts congruently at 2790±30° under a nitrogen pressure somewhat less than 1 
atm [9|. On the other hand, using the above calculated elastic constants, we can derive the 
value of transverse and longitudinal elastic wave velocities, i.e., Vt and vi, and then we can 
get the average wave velocity v m by Eq. (JHJ). From the relation between 6d and v m [Eq. 
(J2J)], we can finally get the value of 9d- Figure [3] (b) shows the pressure dependence of the 
elastic wave velocities. One can see that all the three wave velocities increase with aug- 
menting pressure. Through the study of evolution of Debye temperature, we can calculate 
the melting temperatures using the Lindeman criterion. As is shown in Fig. HI the melting 
points of ThN have been calculated from ambient pressure to 117 GPa. It can be seen that 
the melting temperature of ThN is increased by about 2100 K from ambient pressure to 117 
GPa. 

D. Formation energy 

The formation energy of a specific compound is defined as the difference between the 
total energy of the compound and of its constitutive elements. The composition reaction of 
ThN is as follows, 

Th + -N 2 -> ThN, (6) 
2 
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FIG. 5: (Color online) Calculated phonon dispersion curves (left panel) and corresponding pro- 
jected phonon DOS (right panel) for ThN. 

which yields the following expression for the formation energy: 

£/(ThN) = £(ThN) - (E Th + ^E N2 ). (7) 

In order to calculate the formation energy Ef, the total energy of ThN, a-Th and N2 dimer 
should be calculated. Density-functional theory is known to overestimate the binding energy 
Eb of N 2 dimer, so it will result in an underestimation of the present reaction energy via the 
E^ 2 term. However, this error can be remedied by shifting the energy o 
the experimental binding energy. The experimental of N2 is 9.9 eV 
the Eb is overestimated by about 0.7 eV/N2. The formation energy Ef is —2.14 eV/atom 
with the correction in the En 2 term. This value is somewhat higher than the experimentally 



N 2 so as to give 
3- In the GGA 



meaured enthalpy of formation AfH=— 1.81 eV/atom at room temperature [28]. Since the 
formation energy calculated here does not take into account the temperature, therefore, it 
can only give insight into the stability of the compounds at low temperatures. 

E. Phonon dispersion curve 

In calculating the phonon dispersion curves and the phonon density of states, the 



Hellmann-Feynman theorem and the direct method 29[ are employed. For the BZ integra- 
tion, the 3x3x3 Monkhorst-Pack Appoint mesh is used for the 2x2x2 ThN supercell which 
contains 64 atoms. In order to calculate the Hellmann-Feynman forces, we displace two 
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FIG. 6: (Color online) Dependence of the Helmholtz free energy F{T,V) on crystal volume at 
various temperatures and the locus of the minimum of the free energy for ThN. 

atoms (one Th and one N atoms) from their equilibrium positions and the amplitude of all 
the displacements is 0.03 A. The calculated phonon dispersion curves along T—X — K — T — L 
directions is displayed in Fig. |5j For rock-salt type ThN, there are only two atoms in its 
formula unit, therefore, six phonon modes exist in the dispersion relations. The projected 
phonon DOS for ThN is also plotted in Fig. [51 Because of the fact that the thorium atom is 
heavier than the nitrogen atom, the phonon DOS splits into two parts with an evident gap: 
one part is in the range of 0-4 THz where the vibrations of thorium atoms are dominant; 
the other part is in the domain of 8-11 THz where the vibrations mainly come from nitrogen 
atoms. Our calculated optical frequency at T is about 10.6 THz, and this value is very close 



to the experimental value of 10.3 THz 18|. The phonon dispersion illustrates the stability 



of ThN, and further indicates that our following thermodynamic calculations are reliable. 
F. Thermodynamic properties 

The Helmholtz free energy F in QHA is investigated as follows: 

F(V,T)=E(V) + F ph (V,T) + F ele (V,T), (8) 

where E(V) stands for the ground-state cold energy, F ph (V,T) is the phonon free energy at 
a given unit cell volume V, and F ele is electron excitation energy. Under QHA, the F ph ( V, T) 

12 




300 600 900 1200 1500 



Temperature (K) 



FIG. 7: (Color online) Temperature dependence of the bulk modulus for ThN. 
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FIG. 8: (Color online) Specific heat capacities of ThN. Experimental data from Refs. 3CJ, [3l( and 
theoretical results from Ref. [32] are also displayed for comparison. 



can be calculated from phonon DOS g{u>) by 



POO 

F ph {V ) T) = k B T / g{u)\n 
Jo 



2 sinh 



2k B T 



duj, 



(9) 



where u=u(V) depends on volume and thus Equation ([9]) contains some effect of anhar- 
monics. F e i e in Eq. (8) can be obtained from the energy and entropy contributions, i.e., 
E de — TS e i e . The electronic entropy S e i e is of the form 



S de (V,T) = -k B J n(e,V)[f\nf + (1 - /)ln(l - f)]de, 



(10) 



13 



where n(e) is electronic DOS, and / is the Fermi-Dirac distribution. While the energy E e \ e 
due to the electron excitations takes the following form 

E ele (V, T) = j n(e, V)fede - J** n(e, V)ede, (11) 

where Ef is the Fermi energy. 

The calculated free energy F(V, T) curves of ThN for temperature ranging from up to 
1500 K are shown in Fig. 6. The locus of the equilibrium lattice parameters at different 
temperature T are also presented. The equilibrium volume V(T) and the bulk modulus 
B(T) are obtained by BMEOS fitting. Figure 7 shows the temperature dependence of 
the bulk modulus B. Clearly, the bulk modulus B(T) decreases along with the increase of 
temperature. This kind of temperature effect is very common for compounds and metals. 
Besides, the specific heat at constant volume Cy can be directly calculated through 

Cy = ( W )y = k B I ,(„) (— ) (exp J'I 1)2 ^. (12) 

while the specific heat at constant pressure Cp can be evaluated by the thermodynamic 
relationship Cp—Cy = a v (T)B(T)V(T)T, where the isobaric thermal expansion coefficient 
can be calculated according to the formula ay(T) = y (f^) p - Calculated Cy and Cp of 
ThN are displayed in Fig. El For comparison, the experimental data from Refs 
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31] and 



the theoretical molecular dynamics (MD) results by Adachi et al. 32] are also presented. 
In general, our calculated values of Cy and C p are both higher than the corresponding MD 
results up to 1500 K. However, our calculated Cp, including both lattice and conduction 
electron contributions, agrees well with experimental results in a wide temperature domain 
with a tiny discrepancy of 4 J/mol-K at 1500 K. 



IV. CONCLUSION 



In summary, we have performed systematic first-principles calculations on the structural, 
electronic, mechanical, and thermodynamic properties of ThN. Within the GGA method, 
the ground state structure of ThN can be well produced. Our calculated lattice constant is 
in good accordance with the experimental data, within 0.05% error. Calculated electronic 
density of states show the important contributions of Th Qd and 5/ states to the Fermi- level 
occupation. The Bader effective charges of ThN can be expressed as Th +1 " 82 N~ L82 , which 
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is indicated to be more ionic when compared to UN and PuN. The mechanical analysis has 
been carried out, showing that the rock-salt type ThN is mechanically stable in a wide range 
of pressures. Also, we have presented the melting curve of ThN from ambient pressure to 
117 GPa by utilizing the Lindemann criterion. The calculated phonon dispersion of ThN is 
stable, confirming the dynamic stability. Under the QHA method, our calculated specific 
heat, including both lattice and conduction electron contributions, agrees well with experi- 
mental results. We expect that these calculated results will be useful for the application of 
thorium nitrides in the Generation-IV reactor and nuclear industry. 
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